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NOMENCLATURE 


a  half  width  of  structure  or  radius  of  tank 

a,  \/m2-1 

Aj.Aj  defined  by  eq  9  and  1 1 

b  half  length  of  structure 

b)  1/In  (n+d\ ) 

BUB 2  defined  by  eq  9  and  1 1 

C  volumetric  specific  heat 

d  insulation  thickness 

D  phase  change  thickness  beneath  pipe  insulation 

f  steady-state  solution  geometric  function 

f  steady-state  solution  geometric  function  modified  for  insulation 

f0  value  of  f  on  the  phase  change  interface 

g  steady-state  solution  geometric  function 

g0  value  of  g  on  the  phase  change  interface 

9\  P+1 

h  depth  to  center  of  arbitrary  isotherm 

h0  depth  to  center  of  buried  pipe 

H  h/ri 

H0  dimensionless  depth  to  center  of  phase  change  isotherm  circle 

*  thermal  conductivity 

*12  *  i  /*2 

L  volumetric  latent  heat  of  fusion 

Le  effective  latent  heat  =  L  (I  +k12  /JST  *-ST/ 2) 

m  outward  normal  from  phase  change  interface 

m  m/a 

n  b/a 

u  0(1+b,a)/(l+/3) 

p  distance  along  outward  normal 

P  pla 

p,  /3(rr+4a)/(1+/3) 


r0 

R 

«o 

s 

l 

5 

5 

5 

S, 


heat  flux  from  bottom  of  insulated  structure 

heat  transfer  rate  per  unit  length  of  pipe 

radius  of  arbitrary  isotherm 

(*-*')2  +  (y-/)2+/2 

outer  radius  of  insulation 

radius  of  pipe 

rlr  i 

dimensionless  radius  of  phase  change  isotherm  circle 

distance  along  phase  change  interface 

slo 

area  of  phase  change  interface 

Sla2 

surface  area  with  disturbed  temperature 

ci(Vrf>  „  r 

- - - ,  Stefan  number 


IV 


t  time 

T  temperature 

Tf  freezing  temperature 

Tq  ground  initial  and  surface  temperature 

T p  temperature  of  bottom  of  structure  or  pipe  su'f.ice 

7p  temperature  of  insulation/ground  interface 

u  dummy  variable 

V  volume  of  region  changing  phase 

V  via 3 

x,y,z  Cartesian  coordinates 

x0>-f'0>-?0  Cartesian  coordinates  of  phase  change  interface 


a 


(S 


insulation  parameter  for  buried  pipe 
insulation  parameter  for  other  geometries 


t  (Trro) 

2'  (Tp-T, 


y 

fo 

T'x 

7. 

f 


$0 

v 

v0 

K 

P 


$0 


«* 

«c- 

$B 

?8~ 

«P- 

P 


value  of  i|0  at  center  of  semi-infinite  strip 
value  of  7  for  uninsulated  strip  (a  =  0) 
value  of  7  at  any  location  f 
maximum  or  limiting  value  of  7 


*lr  j 

buried  pipe 

x/a 

other  geometries 

jV'i 

buried  pipe 

■  V® 

other  geometries 

y/a 

V® 

thermal  diffusivily 

Vri 

\',r' 

buried  pipe 

z/a 

other  geometries 

Vri 

buried  pipe 

V® 

other  geometries 

value  of  beneath  center  of  circular  tank 
limiting  value  of  £c 
value  of  $0  at  center  of  rectangle 
limiting  value  of 

value  of  |0  directly  beneath  buried  pipe 
limiting  value  of 

dimensionless  radial  distance  from  center  of  phase  change  isotherm 


v 


fe)  _  dimensionless  time  for  buried  pipe 

r?  Z. 

I 

r 

2<C'5t/_  ,  dimensionless  time  for  other  geometries 

n  a 2 

t*  Vr 

0  polar  coordinate 

T-T0 


0Q 


0, 


Tj-h 

Vrf 

Ty-Tf  f-fp 
Tp-Tf  'l-'o 


T2~To  _  f_ 
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1  thawed  region 
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f  fusion  or  frozen 
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CONDUCTION  PHASE  CHANGE 
BENEATH  INSULATED  HEATED 
OR  COOLED  STRUCTURES 


Virgil  J.  Lunardini 


INTRODUCTION 

Large-scale  exploration  and  development  of  the  northern  regions  of  the  Northern  Hemisphere 
have  stimulated  interest  in  a  number  of  thermal  problems.  Not  the  least  of  these  Is  the  effect  of 
heated  structures  on  underlying  or  surrounding  permafrost.  This  involves  the  study  of  conduction 
heat  transfer  in  media  which  can  undergo  freezing  or  thawing.  Lachenbruch  (1957a,b,  1959)  and 
)umikis  (1978)  applied  linear  conduction  theory  to  the  effect  of  heating  on  permafrost.  Linear 
theory  was  possible  since  no  phase  change  was  considered  despite  the  direct  reference  to  perma¬ 
frost  problems.  However,  if  phase  change  is  introduced,  the  conduction  problem  becomes  non¬ 
linear  and  only  a  few  exact  solutions  exist  for  the  simplest  geometries  and  boundary  conditions 
(Lunardini  1981a). 

Geometries  of  practical  interest,  such  as  those  considered  in  this  report,  do  not  allow  exact 
solutions  of  the  thermal  problem  to  be  found.  Thus  closed  form  solutions,  as  opposed  to  numeri¬ 
cal  evaluations,  have  relied  on  approximate  methods.  In  general  two  approximate  methods  have 
been  quite  fruitful  for  phase  change  conduction  problems:  the  heat  balance  integral  method,  and 
the  quasi-steady  method.  The  heat  balance  integral  method  gives  excellent  accuracy  and  has  been 
applied  successfully  to  simpler  geometries  such  as  semi-infinite  plane  systems  (Goodman  1958, 
Lunardini  1981c,  Lunardini  and  Varotta  1981,  Yuen  1980),  and  pipes  buried  at  infinite  depth 
(Lunardini  1980,  Bell  1978). 

The  quasi-steady  method  is  not  as  rigorous  as  the  heat  balance  integral  method  but  it  can  be 
applied  to  a  wide  variety  of  geometries.  Applications  have  included  uninsulated  buried  pipes 
(Hwang  1977,  Thornton  1976,  Porkhayev  1963),  insulated  buried  pipes  (Lunardini  1981a, b,d, 
Seshadri  and  Krishnayya  1980),  and  three-dimensional  structures  (Porkhayev  1970).  Probably 
the  most  widely  used  calculated  results  are  those  of  Porkhayev  (1970).  This  report  derives  new 
quasi-steady  relations  which  are  applied  to  insulated  geometries  including  semi-infinite  strips 
(roads,  sewers),  rectangular  buildings,  circular  storage  tanks,  and  buried  pipes.  The  quantitative 
results  for  insulated  systems  are  superior  to  those  of  Porkhayev  (1970). 


QUASI-STEADY  METHOD 


L 


The  quasi-steady  approximation  assumes  that  the  temperature  field  varies  successively  from  one 
steady  state  to  another.  Let  us  examine  the  limitations  of  this  approximation.  Consider  an  infinite 


TP 


Figure  1.  Insulated  semi-infinite  strip. 


strip  as  shown  in  Figure  1.  Initially  the  temperature  is  uniform  at  T0  and  the  insulated  surface  of 
the  strip  jumps  to  7"p  at  time  zero.  The  temperature  beneath  the  strip  then  starts  to  vary  with  time. 
The  properties  of  the  material  differ  for  regions  below  and  above  the  lusion  temperature  Tt.  The 
equations  for  the  conduction  heat  flow  problem  are 


a2  t  ,  a2^  ,  a  r, 

3x2  3<r2  Kt 

(11 

T,  (x,0,f)  =Tp  x  *  a  t  >  0 

(la) 

I ,  (xQ,/0,t)  =  Tf 

(1b) 

a  2r2  a2r2  ,  ar2 

3x2  3/2  K  2 

(2) 

T2  (x,0,f)  =  Tq  x  <  -a  or  x  >  +a 

(2a) 

T2  (*  pZ ,0)  -  7 o 

(2b) 

T2  (zr 0 >^0 ■  ^ *  = 

(2c) 

lim  T2  (x,z,t)  -*  T0 
x,z-°° 

(2d) 

If  thawing  is  considered,  the  energy  balance  over  the  surface  S  separating  the  frozen  and  thawed 
zones  is 


'  k7 


(3) 


whece  m  is  the  outward  normal  to  the  interface  and  V  the  volume  of  the  material  changing  phase. 
Equations  1-3  can  be  nondimensionalized  using  the  following  quantities: 
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(5a) 
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1 
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Equations  4-6  cannot  be  solved  exactly,  but  the  diffusion  equations,  eq  4  and  5,  reduce  to  the 
steady-state  case  if  the  Stefan  number  ST  is  small.  In  this  case  it  is  not  necessary  to  solve  the  tram 
sient  conduction  equation  but  only  the  much  simpler  steady-state  conduction  equation.  For  sim¬ 
plicity  assume  that  the  properties  of  the  frozen  and  thawed  media  are  identical.  Then  the  quasi¬ 
steady  system  to  be  solved  is 

d20  ,  d20  -  0 
3f2  3t2 
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(7) 
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(i  -i  ^  f ^  i 

0  (f,  o)  -  |  (7a) 

1 0  f  >  I  or  f  <  -I 

I im  o  (f,  £)  =  0  .  (7b) 


Notice  that  the  initial  condition  cannot  be  satisfied  and  that  the  energy  boundary  equation  is  solved 
by  merely  substituting  in  the  solution  for  0.  It  has  been  shown  (Lunardini  1981a)  that  eq  7-7b  are 
the  zeroth  system  of  a  perturbation  solution  of  eq  4-5.  The  quasi-steady  method  then  consists  ot 
solving  the  steady -state  form  of  the  problem  and  locating  the  phase  change  boundary  with  the  use 
ot  an  equation  similar  to  eq  6.  Clearly  the  accuracy  of  the  quasi-steady  method  depends  upon  the 
magnitude  ot  the  Stefan  number ,  which  is  the  ratio  of  die  sensible  to  the  latent  heat  tor  the  material. 
For  systems  with  a  large  latent  heat  relative  to  the  sensible  heat  it  can  be  expected  that  the  quasi¬ 
steady  approximation  will  be  reasonably  good.  This  covers  soil  sy  stems  and  many  phase  change 
materials  used  tor  latent  heat  storage. 

The  utilitv  ot  the  quasi-steady  methods  lies  in  the  fact  that  for  many  phase  change  problems  the 
equivalent  steady-state  solution  can  be  written  down  immediately  or  easily  tound. 

The  solution  to  eq  7— 7 b,  which  is  well  known  and  will  be  useful,  is 


GENERAL  QUASI-STEADY  RELATIONS 

General  equations  can  be  derived  which  will  be  valid  for  a  class  ol  important  phase  change  prob¬ 
lems.  Assume  that  the  solution  to  the  quasi-steady  form  of  eq  4  is 

0,  =  A,  (r)  +fi,  (r)  /'(?,)■)  .  (9) 

This  will  satisfy  the  quasi-steady  form  of  eq  4  and,  when  combined  with  eq  4a  and  4b,  gives 

t-t0 

01  -  T7T  '  (l0) 

1  fo 

The  function  f  is  the  solution  to  the  equivalent  steady-state  equations  and  f0  is  the  value  of  the 
function  on  the  phase  change  surface  where  the  temperature  is  at  the  fusion  value. 

Also  let 


02=A2(r)+fi2(r)/(r,?).  (11) 

Then  using  by  eq  5a  and  5c, 

02  =  r  02) 

'o 

Once  T0  has  been  determined,  the  total  solution  will  be  known.  Equation  6  will  be  evaluated  at 
only  one  location  on  the  phase  change  boundary: 

(ui 

Tp-Tf  dm  dm  n  dr 
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Figure  2.  Geometry  for  thaw  at  ($0,  {g/ 


The  relationship  between  dp' ,  the  small  amount  ot  material  thawed  during  dr,  and  c/£0  is  shown  m 
Figure  2: 


With  these  relations  eq  1  3  is 

( rP~rf  i  ,  4A  _  2 

Vp-  rf  V1  V  W?o.«o  _  *  dT 


(16) 


This  equation  will  give  the  vertical  phase  change  depth  as  a  function  of  time  for  any  location  f0- 
It  will  be  convenient  to  later  evaluate  this  equation  along  the  plane  of  symmetry,  f  =  0. 

It  is  necessary  to  consider  the  effect  of  the  insulation  since  the  temperature  of  the  ground  sur¬ 
face  T p  is  an  unknown  function  of  time.  There  are  several  ways  tc  handle  the  insulation. 


Method  1 

The  heat  flow  through  the  insulation  will  be  equated  to  the  heat  flow  through  the  thawed  soil 
at  $  =  0.  This  concept  has  been  used  by  Scshradi  and  Krishnayya  (1980)  and  by  Lunardini  (1981b) 


(0,0)  (r„-r„) 

a  1  P  f)  6$  '  d 


(17) 


From  eq  10,  this  yields 


T,-h 


1 _ 

2o  a/  (o,o)' 

i-f0  a* 


(18) 
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whore  a  is  the  insulation  parameter  given  in  the  Nomenclature.  The  general  interlace  equation,  eq 
lb,  iv  now 


' _ I _ , _  +  0\0l\ 

/0-.+2a^S  tr 


2  di-Q 
ir  dt 


The  tinal,  steads -state,  or  limiting  solution  occurs  when  d^Qldr  0  in  eq  19: 


'0"-M2q  ~3? 


ai(Q,Q)  +  /0I  =  0 


a/ (o,o) 


The  heat  tlu\,  into  the  thawed  soil,  at  the  center  of  the  heated  surface  is 


This  can  be  written 


-A’t  (^p-^t)  3/  (0,0) 


a  I  -fg - « 


3T  (0,0)1  aj 


To  evaluate  eq  19  it  is  only  necessary  to  find  the  appropriate,  steady -state,  geometric  lunciion  t. 


Method  2 

Method  2  accounts  for  the  effect  of  the  insulation  by  considering  an  excess  layer  of  soil  with  a 
thermal  resistance  equal  to  that  of  the  actual  insulation.  The  excess  soil  layer  is  only  applied  to  the 
temperature  relations  of  the  thawed  /one.  The  concept  was  introduced  by  Porkhayev  (1963,  1970). 
The  temperature  equations  arc 

T,-T,  I' -I  a 

,2J) 

pa 

'r  'o  9o 

where  f  is  the  steady-state  solution,  with  the  excess  thawed  soil  for  insulation,  and  ?  is  the  usual 
steady-slate  solution.  Equation  3  then  becomes 


s  o  •  *  o 


2  J/df Y  +  /3 A2  dJo 

'  n  »\9f  /  +  \3| /  dt 


The  center  heat  flux  to  the  thawed  soil  is 


6 


)  ( Vrf)  3/-' (0,0) 

0(1  -f'0)  a? 


While  both  methods  approximate  the  steady-state  solution  when  an  insulation  layer  is  present,  it 
will  be  shown  that  the  first  method  gives  more  accurate  results.  A  more  accurate  approach  would 
exactly  solve  the  steady -state  problem,  probably  by  using  complex  variables  with  a  Schwarz-Chris- 
toffel  transformation,  but  such  a  method  has  not  yet  yielded  convenient  solutions. 

The  same  relations  apply  to  the  freezing  case  if  r  =  2fc,  (Tt-Tp)  t/{-o  L)  and  region  1  is  the  fro¬ 
zen  layer. 

It  is  now  possible  to  examine,  quantitatively,  several  practical  geometries. 


SEMI-INFINITE  STRIP 


A  semi-infinite  strip  can  represent  roads,  shallow  rivers,  or  very  long  rectangular  buildings.  The 
geometry  of  the  system  is  shown  in  Figure  1 . 

The  steady-state  solution  is  given  by  Lunardini  ( 1 981  a)— see  eq  8-and  the  geometric  function  f  is 


Method  1 

Using  eq  28  and  eq  1 9  the  phase  change  interface  equation  is 


dr  : 


1^‘tS-') 

+  4  rto 

(tl-ll-) 

1  ,  P' 

[v-4f  f°] 

Thus 


[ ($~2+t/2-  I)2  +  4u2]du 


(?2_£/2_1) 


_ I—  +  1 

fn-1-  ~ 


where  =  -  tan'1  /  — - - - \ 

°  n  \?2+U2-l/ 

The  limiting  interface,  at  t  =  °°,  is 
$0-  =  cotPl  +  >/cot2  p,  -  (?2  -  1) 


(29) 


(30) 


(31) 


where Pj  =  ~j  (rr+4a). 

Since  eq  30  is  valid  for  any  location  f,  it  will  be  convenient  to  evaluate  the  equation  at  f  =  0. 
This  is  along  the  axis  of  symmetry  for  the  geometry.  Then 


r 


/ 


(u2  +  1 )  du 
1 


+ 


!_ 

f(u) 


(32) 
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where  /  =  (2/rr)  cot'1  u,  and  7  is  the  value  ol  along  the  axis  ot  symmetry. 
Equation  31  reduces  to 

>-  =  co1  |(T+T  (?  +  2a,l  ■ 


133) 


Since  the  value  of  tQ  is  constant,  the  value  ol  the  thaw  depth  at  any  f,  7> ,  can  be  found  as  a 
function  of  the  centerline  value  7: 


•<  ^  •  -  - 


(34) 


When  7  <  I,  the  phase  change  at  the  edge  of  the  strip  (f  1.0)  is  zero.  Iheretore  eq  23  lot  the 
surface  heat  flux  is 


2*,  (Vr,> 


4a 

I  +  —  - 

7! 


(35) 


Method  2 

The  equations  for  method  2  will  only  be  written  for  the  axis  of  symmetry .  1  hen  eq  2b  is 
Cy  du 


(3b) 


(l-/')[(r/  +  2a)2  +  l|  «/0(</2+ I) 


where 


^  tan'1  (his  assumes  that  a  thaw  sod  layer  2a  will  compensate  for  fbe  in 

sulation  layer 

^  =  \  Un"'  G)- 

The  limiting  solution  is 

I  _  _ £_ 


(l-/'o)l(T«  +2a)2  +  1 1  <,t0(7£+>) 


(37) 


The  solutions  of  methods  1  and  2  are  identical  for  uninsulated  strips  (a  =  0)  but  they  differ  if 
the  strip  is  insulated.  Since  an  exact  solution  to  this  problem  is  not  available,  the  stead-state  solu¬ 
tions  were  compared  to  a  numerical  finite-difference  calculation.  Figure  3  shows  the  steady-state 
centerline  temperatures  for  a  =  0.6,  0  =  0.7.  Clearly  the  temperatures  predicted  by  method  I  are 
superior  to  those  of  method  2.  Note  also  that  the  depth  of  the  thaw  zone  predicted  by  method  1 
is  considerably  more  accurate  than  that  of  method  2.  All  of  the  quantitative  results  were  calculated 
with  method  1  due  to  its  superior  accuracy. 

Equation  32  can  be  evaluated  exactly  if  p  =  0.  For  this  case 

T=(’  +  t)(¥  +y)corly+  f  mi+72)^  yj  <38> 

Equation  32  was  evaluated  numerically  and  plotted  as  Figures  A1  -A1 2.  Equation  3 1  is  plotted 
as  Figures  A13-A15  for  various  values  of  f. 

The  temperatures  in  the  thawed  and  frozen  zones  can  be  found  with  eq  10  and  12. 


? 


Figure  3.  Steady-state  temperature  comparisons  be¬ 
neath  center  of  strip. 


RECTANGULAR  BUILDING 


While  the  semi-infinite  strip  can  represent  a  building  with  a  very  large  length/width  ratio,  it  is 
useful  to  have  solutions  for  any  aspect  ratio.  The  geometry  for  such  solutions  is  as  shown  in  Fig¬ 
ure  4. 

The  steady-state  solution  is  given  by  Lunardini  ( 1981a)  and  the  geometric  factor  is 


f  = 


2rr 


-i  (ill!  (v*n) 

tV^+(f  +  1)2  +  (17  +  17)2 


-  tan" 


(f-D(r?  +  n) 


«  \/$2  MM)2  Mb+n)2 


-  tan 


-I 


(£  +  1)  iv-n) 

i  nA2  +  (r + 1  )2  +(t?-n)2 


+  tan"' 


_ (?-  1)  (v-n) 

i;\/S2  +(f  -  D2  +  (b-n)2 


(39) 


where  n  =  bja  is  the  aspect  ratio. 

The  calculations  will  be  carried  out  along  the  axis  of  symmetry,  J  =  rj  =  0,  and  the  geometric 
function  can  be  written  as 


Equation  19,  for  the  phase  change  interface,  is  now 
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Figure  4.  Geometry  for  a  rectangular 
building  on  the  ground  surface. 


| «2  (u-  +  1  +  /;2)  +  /;2|  \Ju~  +  1  +/72  du 


(2  u2  +  1  +  ri2 


where 


(41) 


As  the  aspect  ration  becomes  large,  eq  41  will  reduce  to  eq  32  for  the  infinite  strip.  Thus  the  in¬ 
finite  strip  is  represented  by  eq  41  when  the  aspect  ratio  is  large. 

The  limiting  or  steady-state  solution  is 


Numerical  quadrature  of  eq  41  leads  to  the  plots  given  by  Figures  A16-A51  for  n  -  1,2,3.  The 
limiting  steady-state  values  arc  given  in  Figures  A52-A60. 


CIRCULAR  TANK 

Storage  tanks  are  often  used  in  cold  climates.  The  solution  for  a  circular  tank,  shown  in  Figure 
5,  follows  in  the  same  way  as  for  the  other  geometries. 

The  transient  solution  for  the  temperature  in  a  semi-infinite  medium  initially  at  7"q  after  a  sur¬ 
face  area  S,  is  disturbed  with  a  temperature  Tp  is  given  by  Lachenbruch  (1957a)  as 


where  =  (x-x')2  +  (y-y ’)2  +  /2  and  x',y'  =  coordinates  of  dA  in  . 
The  steady-state  solution,  /-«■>,  reduces  to 
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Figure  5.  Circular  tank  on  ground 
surtac  c. 


Equation  44  can  be  integrated  for  certain  simple  geometries  such  as  the  inlinite  strip  or  the  rectang¬ 
ular  area  already  discussed.  A  general  solution  for  the  circle  is  not  available  but  the  temperature 
along  the  /-axis,  x  =  y  =  0,  can  be  written  immediately: 


T-T0 

Vr0 


=  1  - 


z 


(45) 


The  geometric  (unction  for  this  case  is  then 


m  -- 1  - 


(4b) 


Equation  19  for  the  phase  change  interface  depth  can  then  be  written  for  the  area  along  the  axis 
of  symmetry 


T  =  - 


2 

rr 


(u2  +  1 ) 3/2  du 

I _ L _ +  1 _ 

0  f{u)~  l-2o  f(u) 


where  t\u)  =  1  -  —  -  — 
Vu2  +  1 


(47) 


(48) 


Equation  47  is  plotted  in  Figures  A61-A72  for  various  values  of  a  and  ff.  The  limiting  solution, 
plotted  as  Figure  A73,  is 


oo 


(49) 


11 


The  single  phase  solution  (0  -  0)  can  be  written  explicitly  as 


|4  _ _  _ 

T  z  1 Y  +  ^  +  f  kc  (2^c  +  5)  v  S2  +  I  +  3  In  (£t  +  +  I ) 


(50) 


BURIED  PIPE 


Buried  pipes  are  frequently  used  to  convey  mass  and  energy  in  cold  regions.  The  geometry  of 
the  insulated  buried  pipe  system  is  shown  in  figure  6.  The  steady -state  solution  for  uninsulated 
pipes  is  given  by  Carslaw  and  laeger  (1959)  and  for  insulated  pipes  by  Lunardini  (1981b,  d).  The 
geometric  function  is 


as. « 


f2  +  U  +  «il2 


t2  +  u -a,  i 2 


(51) 


The  temperature  functions  are  given  by 
f-fn 


T,-T{ 


Tp-Ti  I  +b|  a-^0 


(52) 


02 


Tf  -  T  n 


'f-'O  '0 

The  general  equation  for  the  radius  of  the  region  that  has  changed  phase  is 

dRn 


rr/2 

r 

30, 

302 

/ 

dp 

V 

-rr/2 

dO  -  -n 


J  p-Rr 


dt 


(55) 


(54) 


From  eq  52  and  53  this  is 


Figure  6.  Geome  try  of  insulated  buried  pipe. 
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rr/2 

/ — ! -  r  /» 

' 0 /  7^2  \3p/p-Ro 


If  the  heat  flux  is  assumed  to  be  constant  around  the  entire  phase  change  interface,  using  the 
value  at  the  bottom  of  the  pipe  where  0  =  rr/2,  then 


(bf\  a\b\ _ 

W.,Rn  "  /?o(«o  +  /?o 


From  eq  55  the  phase  change  radius  is  governed  by 
u ’  {pj  +  ~\)du 

a\b\T~-  J  ~  1  x 


(p2  + 

I  " 

I  ,  Of- /0  /q 


f'o  3  b\  ln  (pi  +  \'p\-  I  ) 
uz 

Equation  56  can  also  be  written  as 

sD 


(^4  -<7j)  du 


r  ...  i  *  1 

f!'+<,l] 

\u-(ly) 

In 

(~J 

The  solution  to  cq  57,  when  0  =  0,  is 

/  F  (  n  4-n  . 


(t  * £)  *«l?  M  (i  ■  »t)I 

-  y  &  -*?) +  ^ In fe)  -  H  ln  (Jbr) 


where  g  |  =  pi  +  1. 

The  limiting  depth  can  be  evaluated  from 

(/t+<7,)M  +  1 

*P”  '  (ai+o,)m-1 


where  M  =  (1  +  t  or) 


In  this  equation  denotes  the  depth  to  the  bottom  of  the  thaw  interface  on  the  plane  of  sym¬ 
metry  where  J0  =  0.  Solutions  to  eq  57  are  plotted  in  Figures  A74-A91  where  D  is  the  thickness 
of  the  phase  change  beneath  the  center  of  the  pipe. 
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The  heat  How  rate  per  unit  length  ot  pipe  to  the  surface  is  given,  at  any  time,  by 


2tt k]b]  (Tp  -  Tt) 
I  +6,  a-  /0 


(60) 


where 


For  hot  oil  pipes  buried  in  permafrost,  the  sensible  heat  is  considerably  larger  than  the  latent 
heat.  Therefore  the  Stefan  number  is  quite  large  and,  since  the  quasi-steady  method  is  based  on 
systems  where  the  Stefan  number  is  near  zero,  these  cases  arc  a  severe  test  of  the  quasi-steady  equa¬ 
tions.  In  order  to  account  for  the  sensible  heat,  an  effective  latent  heat  parameter  can  be  used, 
defined  as 


l.JL~  1  +Ci,A,,>(351  +St/2.  (61) 

CONCLUSIONS 

The  problem  ot  conduction  heat  transfer  with  phase  change  has  been  solved  approximately  using 
the  quasi-steady  method.  The  solutions  for  the  buried  pipe  have  been  compared  to  numerical  solu¬ 
tions  and  the  quasi-steady  results  were  found  to  differ  by  no  more  than  ±20%  for  most  cases 
(Lunardini  198 Id).  The  pipe  solutions  also  agreed  very  well  with  a  heat  balance  integral  solution 
at  infinite  burial  depth  (Lunardini  1980).  Since  these  comparisons  were  for  Stefan  numbers  greater 
than  1.6,  it  can  be  expected  that  the  accuracy  in  general  will  be  better  than  ±20%  because  the 
quasi-steady  method  improves  as  the  Stefan  number  decreases. 

The  method  of  calculating  conduction  heat  transfer  with  phase  change  presented  in  this  report 
has  been  shown  to  give  more  reliable  results  than  the  widely  used  graphs  of  Porkhaycv  (1970). 

The  quasi-steady  method  is  extremely  useful  because  it  can  be  applied  quite  easily  to  a  number 
of  practical  cases,  and  the  results  can  be  presented  in  a  compact  form  as  shown  by  the  graphs. 

These  graphs  are  felt  to  be  acceptable  for  engineering  estimates  if  accuracies  on  the  on  the  order  of 
20%  are  adequate.  For  insulated  systems,  where  the  phase  change  is  expected  to  be  more  limited, 
the  graphs  should  be  even  more  accurate. 
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APPENDIX  A.  CALCULATED  DESIGN  GRAPHS. 


FIGURE 

A1-A12  Phase  change  beneath  center  of  infinite  strip. 
A13-A15  Limiting  phase  change  for  infinite  strip. 

A16-A27  Phase  change  beneath  center  of  rectangle,  n  =  1 . 

A28-A39  Phase  change  beneath  center  of  rectangle,  n-  2. 

A40-A51  Phase  change  beneath  center  of  rectangle,  n  =  3. 

A52-A60  Limiting  phase  change  for  rectangles. 

A61-A72  Phase  change  beneath  center  of  circular  storage  tank. 
A73  Limiting  phase  change  for  circular  tank. 

A74-A91  Phase  change  beneath  buried  pipe. 
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